Decay of a quantum vortex: test of non-equilibrium theories for warm Bose-Einstein condensates 
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The decay of a vortex from a non-rotating high temperature Bose-Einstein condensate (BEC) is modeled us- 
ing the stochastic projected Gross-Pitaevskii equation (SPGPE). In order to exploit the tunability of temperature 
in SPGPE theory while maintaining the total atom number constant, we develop a simple and accurate Hartree- 
Fock method to estimate the SPGPE parameters for systems close to thermal equilibrium. We then calculate the 
lifetime of a vortex using three classical field theories that describe vortex decay in different levels of approx- 
imation. The SPGPE theory is shown to give the most complete description of the decay process, predicting 
significantly shorter vortex lifetimes than the alternative theories. Using the SPGPE theory to simulate vortex 
decay for a trapped gas of 5 x *^Rb atoms, we calculate a vortex lifetime t that decreases linearly with 
temperature, falling in the range 20s> t >1.5s corresponding to the temperature range 0.78Tc < T < 0.93Tc. 
The vortex lifetimes calculated provide a lower bound for the lifetime of a persistent current with unit winding 
number in our chosen trap geometry, in the limit of vanishing vortex pinning potential. 



Bose-Einstein condensates present a unique opportunity for 
developing ab initio theory that can be directly compared to 
experiments. However, in the area of non-equilibrium dynam- 
ics most quantitative comparisons have been limited to the 
zero temperature regime, where the Gross-Pitaevskii equation 
(GPE) provides a comprehensive treatment. 

An area of particular interest occurs in the non-equilibrium 
dynamics leading to the nucleation of vortices into a conden- 
sate, either as it is stirred or as it grows out of a rotating vapor, 
and the subsequent formation of a vortex lattice 1 1 1. Two types 
of theoretical model have been used to describe this process: 
(i) Gross-Pitaevskii treatments in which a damping parameter 
is introduced 12,^, (ii) Hamiltonian classical-field methods, 
in which the nonlinearity within the Gross-Pitaevskii equa- 
tion and appropriately chosen random initial conditions al- 
lows a stirred condensate to relax to a vortex lattice ||4}{6|. 
While this theoretical work (note also the recent works us- 
ing other approaches |T, 'HI) has established qualitative un- 
derstanding of the processes that lead to vortex lattice forma- 
tion, the timescales predicted for lattice formation are qual- 
itatively incomparable with experiments. A general theory 
of vortex formation is further complicated by the variety of 
nucleation scenarios (a number of distinct mechanisms have 
recently been observed [9j J_()J), and the challenges in devel- 
oping a fully time dependent non-equilibrium theory of the 
interacting Bose gas that can account for all of the relevant 
physics at finite temperature. Indeed, the frontiers of many 
body physics ifTTl . and the application of generalized classi- 
cal field theories fT2l to non-equilibrium fTJl, finite tempera- 
ture fT4\, and critical phenomena 1 15 1 in systems of ultra-cold 
Bosons are fields of much current interest. 

In this paper we examine a well-defined process suitable 
for testing finite temperature theories: the decay of a singly 
quantized vortex in a stationary Bose-Einstein condensate at 
finite temperature. From general arguments it is clear that 
such vortices are metastable in the absence of external rota- 
tion lfT6l ITtI . and will decay in the presence of dissipation. 
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FIG. 1. Schematic of different vortex decay processes in two dimen- 
sions, (a) Under Hamiltonian evolution, an eigenstate of the GPE 
containing a displaced vortex will evolve so that the vortex follows 
trajectories of constant radius about the centre of the condensate. At 
finite temperatures, the presence of a non-rotating thermal cloud cre- 
ates a damping mechanism for vortex motion, (b) Dissipative decay: 
in the absence of noise, a damped, off center vortex will decay to the 
edge of the condensate along a regular spiral trajectory, (c) Diffu- 
sive decay: thermal fluctuations lead to stochastic vortex dynamics 
resembling Brownian motion. 



Previous work on this problem has included: a variational 
treatment based on the formalism of Stoof and coworkers 1 18], 
a pure classical-field treatment |19|, simulations using the 
Zaremba-Nikuni-Griffin formalism |20|, and using a damped 
GPE ETl . We also refer to Ref. |22| for a self consistent 
meanfield treatment, which has been limited to low temper- 
atures, with the prediction that the vortex precesses but does 
not decay. 

In the present work we use the stochastic projected Gross- 
Pitaevskii equation (SPGPE) |'7] |23]| to study vortex decay. 
We show, by neglecting various processes, that this formalism 
can be reduced to the damped GPE and classical field treat- 
ments |fT9ll2TI . In more detail, our approach is a c-field the- 
ory (see IIT2I ) in which the modes of the system are divided 
into two regions according to occupation: appreciably occu- 
pied (low energy) modes constitute the coherent region and 
are described by a classical field theory based on the Wigner 
phase space representation; the remaining sparsely occupied 
modes constitute the incoherent region which behaves like a 
thermalized reservoir and is treated semi-classically. In the 
full SPGPE treatment the dynamics of the coherent region are 
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simulated according to a stochastic differential equation that 
includes damping and noise arising from the incoherent re- 
gion [i?]. By neglecting this noise a damped projected Gross- 
Pitaevskii equation (DPGPE) treatment is obtained (similar to 
that in [20, ,21J). By neglecting both the noise and damping 
from the incoherent region (i.e. all coupUng to the incoherent 
region) the projected Gross-Pitaevskii equation (PGPE) clas- 
sical field treatment is obtained (similar to that in 1 19|). 

By comparing simulations of the SPGPE, DPGPE, and 
PGPE, we are able to understand the relationships between 
these theories and to investigate various physical processes 
that contribute to vortex decay. In the absence of any dissi- 
pative/thermal effects the vortex can only precess in the sys- 
tem and it is unable to leave (Fig. [TJa)). The damping ef- 
fect of thermal cloud leads to a dissipative decay of the vor- 
tex (Fig. [TJb)). However, the associated noise from the ther- 
mal cloud interaction will also cause diffusive dynamics of 
the vortex (Fig. [ijc)). Within the c-field formalism it is not 
correct to solely attribute the dissipative and noisy features 
of condensate evolution to the incoherent region, indeed the 
coherent region itself is a multi-mode system with its own 
intrinsic dissipative and fluctuation characteristics. Thus the 
relationship between the important processes and the various 
theories is not immediately clear. The results of these distinct 
formalisms reveal that the SPGPE predicts the shortest vortex 
lifetime by a significant factor, suggesting that the other ap- 
proaches neglect important physics and are not useful quanti- 
tative models. 

Having established the importance of the method, we ap- 
ply the full SPGPE to the study of vortex decay for an experi- 
mentally realistic system over the temperature range 0.78Tc < 
T < 0.93Tc, for traps with spherical and oblate geometries. 
We present a technique for effectively estimating the SPGPE 
simulation parameters to obtain initial (equilibrium) states 
with a fixed total number as the temperature is varied. We de- 
fine a vortex lifetime from the vortex decay trajectories, and 
characterize its variation with increasing temperature, finding 
a hnear decline in the regime studied. A comparison of vor- 
tex lifetimes for two geometries suggests that vortex bending 
modes associated with Kelvin waves provides an important 
mechanism of vortex decay that is controllable through the 
degree of trap oblateness. 

The angular momentum dynamics indicate regimes of both 
dissipative and diffusive decay of the vortex, where the latter 
regime is a consequence of thermal fluctuations from classi- 
cal and quantum modes. The dissipative and diffusive vortex 
decay mechanisms are illustrated schematically in Fig[T] 

I. STOCHASTIC PGPE THEORY AND MODELS 

The derivation of the SPGPE for high temperature Bose 
gases has been given elsewhere ||23]| . Here we briefly sum- 
marize our approach ifTOllTSl . A dilute Bose gas in the cold- 
collision regime is described by the second-quantized many 




FIG. 2. (a) Separation of the high temperature Bose gas using the 
single particle energy cutoff tcut, into coherent (C, with energies 
< Ecut) and incoherent (I, with energies Sli > ecut) regions, (b) 
The I - C collisions that lead to Bose-enhanced growth or loss: dur- 
ing a coOision of two I atoms, one atom absorbs most of the collision 
energy, allowing a stimulated transfer of the other into the coherent 
region, C. The time-reversed process also occurs. 



body Hamiltonian 

H = j d?^ i^^(x)i7,pV(x) 

+ ^ yd3x^t(x)V,t(x)V'(x)V'(x), (1) 

where the single-particle Hamiltonian and trap potential are 

i^.p = -^ + nx), (2) 
nx)=|K(.^ + y^)+c..V] (3) 

and u = Airfi^a/m, with m the atomic mass and a the s-wave 
scattering length. 

As shown schematically in Fig.|2] the states of the trapped 
system are divided into what we refer to as the coherent re- 
gion (C) of states with energy beneath the cut-off Ccut, and 
the remaining incoherent region (I) of high energy, quasi- 
equilibrium states. The C modes are assumed to be, at least 
moderately occupied and hence well-approximated by a clas- 
sical field. This identification is key to the utility of the vari- 
ous c-field techniques as it allows the full coherent dynamics 
of the low energy portion of the system to be treated non- 
pert urbatively. 

The field operator is decomposed as 

V'(x) = ^c(x) + ^i(x) (4) 

where -(AcCx) = 7'c{V'(x)} and ■0i(x) = '(/'(x) — V'c(x), are 
the field operators for the coherent and incoherent regions, re- 
spectively. Our precise definition of the regions is made us- 
ing a projection operator, Vc, defined by its action on a wave 
function /(x) according to 

Vc {/(x)} = / -^"(^O/Cx'), (5) 

where the 0„(x) are eigenstates of Hgp, i.e. e„0„(x) = 
Hsp4>n{x), and the restriction to the C region in the summa- 
tion is given by {n : e„ < Ccut}- 
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A. Semi-classical treatment of the I region 

In the SPGPE treatment the I region modes act as a thermal 
reservoir for the low energy coherent region and are treated 
using master equation techniques. The situation we consider 
here, of a vortex damping from a partially condensed gas, is 
near equilibrium and it is adequate to assume that the I region 
(which thermalizes rapidly compared to the C region) is in 
equilibrium with a well-defined temperature T and chemical 
potential /i, and is well-described by a single particle Wigner 
function 



Fi(x,k) 



exp [(/iu;(x, k) — fi)/kBT] — 1 ' 



where 



2u2 



huj(jx., k) = 



2m 



(6) 



(7) 



By virtue of our cutoff this description is restricted to the I re- 
gion, and any quantities computed from (|6]l should only be 
integrated over the incoherent region of phase space Ai = 
{(x,k) : huj{x,'k) > Ecut}- For example, the number of I 
region atoms is given by 



A^i-<?3(e'^^,/?ecut)(/3to)"', 
where Q = (w^a;^)^/'^, and 

pec ^ 



(8) 



(9) 



1=1 



is the incomplete Bose-Einstein function 

B. SPGPE theory of the C region 

By formally tracing over the I region, we arrive at a de- 
scription for the dynamics of the C region field in which the 
influence of the I region appears as damping (7) and complex 
noise {dW^) terms. The evolution equation has the form 

d^c(x,t) =Pc{ - ^LcM^,t)dt 

(10) 

which we refer to as the simple growth SPGPE. The operator 
Lc is the Hamiltonian evolution operator for the C region 



icV^c(x) = (H.p + w|7^c(x,t)|2) ^c(x), 



(11) 



and the complex noises are given by 

{dW;{x,t)dW^{x\t)) = 2-f6ci^,^')dt, (12) 
{dWj{yi,t)dW^{x',t)) = {dW*{x,t)dW*{yi',t)) ^ 0,(13) 

where (5c (x, x') = J^nec '/'n(x)'/'*i(x') is a delta-function in 
C. This result is systematically derived in Ref. Q) using the 



truncated Wigner approximation and master equation tech- 
niques, with the notable feature that the damping and noise 
parameters are calculable from our description of the I region 
for a system close to thermal equiUbrium, giving the expres- 
sion 



7 = 7o< 



In (1 - e'^^'^-'^-')) + 



„2/3(M-£cut) 



E 

r=l 



(14) 



4:m{akBT)^ /nh^, and ^[x,y,z] is the Lerch 



where 70 
transcendent. 

The second line of the SPGPE ([TOjl is directly responsible 
for condensate growth from scattering between two I region 
atoms, as shown schematically in Fig. |2|b). ^ and T are re- 
spectively the chemical potential and temperature of the ther- 
mal reservoir of particles in the I region, defined in Eq. (|6|. 



By using the spatially invariant form ( 14 1 for the damping, we 
have neglected a weak spatial variation arising at the edge of 
the C region |7 1. 

Additional terms from the full SPGPE theory that represent 
number conserving scattering processes between atoms in C 
and I 1231 have also been omitted. These processes do not 
lead to condensate growth, and are weak in near equilibrium 
situations [12|. The numerical implementation of Eq. ([TOjl is 
relatively straightforward, with individual trajectories being 
only slightly more computationally intensive than integration 
of the PGPE. Our method of numerical stochastic integration 
has been given in a number of places; we refer the reader to 
Refs. ElEllSll for further details. 

We note the related stochastic GPE approach of Stoof and 
coworkers ITS', "241 (also see |251), which differs from our ap- 
proach most significantly through the absence of an explicit 
projector. A comparison of these approaches can be found 
in Ref. 1*261. It is our opinion that the projector is a crucial 
element of the reservoir theory that cannot be simply chosen 
arbitrarily, as is extensively discussed in Ref. |12J. 

The SPGPE description provides a grand canonical de- 
scription of the C region so that, irrespective of the value 
of 7, in equilibrium it samples microstates with probability 
P(V'c) oc exp {—Kq/ kgT), where Kq = Hq — fJ-N^ is the 
grand-canonical Hamiltonian and 

He = y d'x^*(x) (i/,p + ||^c(x)p) Vc(x), (15) 

Nc= f d3x|^c(x)P, (16) 



are the expressions for the c-field energy and number respec- 
tively. We now consider two simplified theories that can be 
obtained from the SPGPE by neglecting various effects. 
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C. Damped PGPE treatment 

A damped PGPE (DPGPE) is obtained by neglecting the 
noise term in the SPGPE ([T0|, and is of the form 

■f^di^c^ =Pc {ic^c(x, t) ~ ^^{Lc - A^)^c(x, t) 

(17) 

While the SPGPE serves to minimize the grand free energy 



of the C region, evolution according to equation ( 17 1 instead 
minimises Kq: 



dKc 

dt 



27 



d3x|7'c{(M-ic)V^c(x)}|^ (18) 



i.e. thermal fluctuations are damped out, and the equilib- 
rium solution is the zero temperature ground state, satisfying 
VcLctpc = iJiipc- An initial state containing a vortex will 
decay to a vortex-free state (provided the I-region is rotating 
sufficiently slowly), on a timescale set by 7. 

To date several numerical studies of vortex nucleation and 
dynamics have been based upon a phenomenological damped 
Gross-Pitaevskii equation ||2l [3] |22l. The description was 
first introduced in fT\ on the grounds that it gives the cor- 
rect equilibrium solution, and that it follows naturally from 
the damped GPE approach of ESl for a cloud in equilibrium 
at rotation frequency VI. This approach received further justi- 
fication via a stochastic treatment of the dissipative interaction 
of a BEC with a rotating thermal cloud ||7]|23|. As shown by 
Kasamatsu et al. |27 1, the equation of motion may also be ob- 
tained from the generalized finite temperature GPE formalism 
of Zaremba, Nikuni, and Griffin (ZNG) 1291. 



D. PGPE treatment 

By setting 7 (in addition to dW^ = 0) we arrive at an 
isolated PGPE description of the C region 



-9V'c(x,t) 

' dt 



=Vc{LcM^.t)}- 



(19) 



The PGPE is also formally obtained as Hamilton's equation 
of motion for the c-field Hamiltonian ( [T5] l: 



dt 



SHc 



(20) 



where the restriction to the C region arises naturally in the 
projected functional differential formalism of Ref. 123)1 . An 
interesting feature of the PGPE theory is the emergence of 
ergodicity from Hamiltonian evolution lfT2l (particularly for 
the finite temperature regime we consider, where the C region 
contains many appreciably occupied modes). While the PGPE 
approach may appear to be a simplification of the DPGPE, it 
provides quite a distinct physical description. It is a micro- 
canonical model of the C region in which both Hq and Nq 
are conserved and, by virtue of 7 being zero, the full comple- 
ment of coherent region excitations (and their intrinsic noisy 
effect on the C region dynamics) are retained. 



Two recent applications of a rotating frame PGPE (7) have 
studied vortex lattice formation [6 |, and symmetry breaking 
associated with the single vortex in equilibrium l30l . Other 
treatments of the isolated coherent region have been per- 
formed, typically without a projector and have adopted the 
classical field moniker Important earlier works are those of 
Schmidt et al. |19|, who considered the decay of a vortex in 
a purely Hamiltonian GPE evolution of a thermal Bose field 
containing a vortex, and of Parker and Adams lH] [3Tj who 
studied the formation of a vortex lattice in a stirred, zero tem- 
perature BEC. 



II. STATE PREPARATION AND VORTEX DETECTION 

In our simulations the same initial conditions are used for 
the SPGPE, DPGPE and PGPE calculations, and are aU pre- 
pared using the SPGPE to obtain an finite temperature equi- 
librium state. A central vortex is then generated in the state, 
and its position is detected as a function of time. In this sec- 
tion we discuss the procedure used to generate these states and 
distinguish the vortex from surface fluctuations. 



A. Efficient sclieme for preparing initial states at fixed A'^t 

In many applications it is desirable to perform simulations 
in which the temperature is changed but the total atom num- 
ber in the system, iVr (and hence critical temperature Tc), re- 
mains constant. Equilibrium states in the SPGPE description 



are generated by evolving Eq. ( lOl with choices for T and /i, 
but as the relationship between ji and iVr is nonlinear, obtain- 
ing the desired iVr (which can only be found a posteriori) can 
require running a large number of simulations to determine 
the appropriate ji value by trial and error Another parame- 
ter that must be determined is the value of ecut used to define 
the C region. The validity of the SPGPE approach requires 
that the highest mode in the C region be appreciably occupied 
(usually with a occupation in the range ricut = 1 — 5) and to 
ensure this value of Ccut will also change with T where neces- 
sary (e.g. see 1*3211). Because the projection is made using the 
single particle basis, we also require a sufficiently large cutoff 
to describe the interacting modes of the system. In practice 
this requirement is satisfied if ecut ^ 2/i. 

Here we develop a simple mean field method that provides 
useful estimates for the choices of the parameters /i and ecut 
for a system of atoms at temperature T. 

In our approach the condensate is treated within in the 
Thomas Fermi approximation, with the number of conden- 
sate atoms (iVo) determining the system chemical potential. 
The thermal atoms are then treated using the semi-classical 
Hartree-Fock approximation, i.e. with energies 



^hf(x,p) = ;i^ + y(x)+2u7io(x) (21) 
Zra 

where no(x) is the condensate (Thomas-Fermi) density. The 
Hartree-Fock density of states can not in general be obtained 
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FIG. 3. Estimated and SPGPE parameters as a function of input tem- 
perature using an input value of A'^t = 1 x lO*" for a spherical system 
with trap frequency a; = 27r x 10 s~^. (a) Hartree-Fock estimates of 
the cutoff energy ecut for a range of cutoff mode occupations ricut- 
The dashed lines give the curves for tcut and fi used in (b) Total 
atom number from SPGPE simulations evolved to equilibrium using 
HF estimated parameters (circles), and fixed jj, (squares), /i is held at 
the value giving TVt at T = 34 nK. The solid line is the target value 
of A^T. 



analytically, but it can easily be handled for the case of a har- 
monically trapped condensate, as shown in Ref. |33|. Fol- 
lowing this approach we have obtained an efficient numerical 
algorithm for estimating /i and ecut at a given temperature and 
total atom number The details of this procedure are given in 
Appendix [a] 

Results from the Hatree-Fock calculation for /i and Ccut for 
a system of 1 x 10^ ^^Rb atoms are shown in Fig.|3](a), reveal- 
ing the typical way in which these parameters need to change 
with temperature to ensure fixed total number We see that the 
value of Ccut for ricut = 2 is sufficiently large that the require- 
ment Ccut is well-satisfied in the temperature regime 
O.TTc < T Tc. It should be noted that the Hartree-Fock 
predictions only depend on the trapping frequencies through 
u}, so that these parameters can be applied to any system of 
equal geometric mean, provided a three dimensional descrip- 
tion is appropriate. 

Because of the simplicity of the Hartree-Fock theory it is 
necessary to test that its predictions are accurate. In what fol- 
lows we present SPGPE results using the Hartree-Fock esti- 
mates for /i and ecut to verify the accuracy at which the desired 
total number is obtained for various temperatures. To do this 
for each set of parameters we determine the incoherent region 
population, A'l, according to Eq. (j8]l, and from the equilibrium 
state obtained from the SPGPE simulation using the value of ji 
and ecut predicted by Hartree-Fock theory, we determine Nq, 



TABLE I. Comparison of input and resulting SPGPE reduced tem- 
peratures for a spherical system with a;r = Stt x 10s~^ and input 
atom number A^t ~ 1 x 10^. 



and thus the total number iVr = A^c + Ni. To obtain equilib- 
rium states, the SPGPE [Eq. ( [TO) i] is evolved with the param- 
eters determined using the Hartree-Fock scheme. Note that 
the grand canonical equilibrium solution is independent of the 
size of 7, allowing us to find equilibria more rapidly using a 
dimensionless damping H-j/KbT ^ 0.3 that is much larger 
than typical dynamical values {Hj/KbT ^ 10^"^). We then 
find that the c-field reaches thermal equilibiium after prop- 
agation for much less than a single radial trap period, when 
initialized from a Thomas-Fermi initial state. 

Our results are calculated for a Bose gas in a spherical trap 
with frequency a; = 27r x 10 s^^ and Nt = 1 x 10^ atoms. 
Figure [5] (b) shows the SPGPE results for A't over a range 
of temperatures and reveals that the HF estimated parameters 
generate equilibrium SPGPE solutions with total atom num- 
ber in close agreement for the desired number For compari- 
son we also show the effect of holding /i constant and increas- 
ing temperature in Fig.[3](b). In this case the total number of 
atoms increases quite rapidly, and hence the critical tempera- 
ture also increases (/cbTc = 0.94toA/i/^). Table [ijcompares 
the two methods, for the results of Fig. [4] The SPGPE results 
using the HF parameter estimates effectively and accurately 
span the temperature range where the SPGPE is valid. In 
comparison, the fixed /x case approaches the transition slowly 
since increases. We have also verified that the method 
works over a range of geometries. The c-field occupation at 
the cutoff is typically slightly smaller than the input value. 
For example, we find that using an input value of ricut = 2 
to the HF-estimate of Ccut results in SPGPE solutions with 
"cut ~ 1-5. 



B. Preparing a vortex state 

Finite temperature states containing a vortex line that lies 
on the cylindrical symmetry axis of the trap (z), are induced 
by imposing a vortex phase pattern on samples of the SPGPE 
equilibrium state, ipc{^) ?Ac(x)e*'^^''\ where 



(x) = arctan 



(22) 



Imposing a phase consistent with a single quanta of circula- 
tion has the effect of only coupling to the coherent region, so 
that the thermal fraction above ecut and the high energy ther- 
mally occupied modes of the c-field region are not changed. 
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(a) {b) (c) 
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FIG. 4. Vortex dynamics using the DPGPE and POPE in a system with a 4:1 geometry. The vortex is most easily observed from the phase 
singularity (circle), and the center of the condensate is marked by a cross (+). (a): DPGPE vortex evolution for T/Tc = 0.78, where the vortex 
has an initial offset from the origin of the condensate. The initial noise is damped out within ~ 1.2s, after which the vortex decay is governed 
by 7. (b): PGPE vortex evolution for T/Tc — 0.78, with an initially centered vortex. Despite initial thermal fluctuations kicking the vortex off 
the center, the size of the thermal component is too small at this temperature to see a finite lifetime, (c): PGPE vortex evolution for a system 
at T/Tc ~ 0.93 where the initially centered vortex has a finite lifetime due to increased classical fluctuations at this temperature. 



Thus the system is initialized as a finite temperature equilib- 
rium state with non-rotating thermal component, containing a 
central, axially aligned, singly charged vortex. 

C. Vortex detection 

To unambiguously detect the a vortex as it migrates to the 
condensate edge, we use an adaptive method that compares 
the vortex signal with surface fluctuations. In detail, the curl 
of the velocity field w(x) = Vxv(x),is used as our basic sig- 
nal. Taking a slice through the c-field at z = 0, a central vor- 
tex generates a significant signal localized at the vortex core, 
which is unambiguously resolvable; we denote the vortex dis- 
placement from the z-axis by However, near the conden- 
sate boundary, the vortex is more difficult to distinguish from 
phase singularities in the evanescent matter field caused by 
the stochastic noise term in the SPGPE. In experiments this 



would correspond to a loss of core visibility in absorption im- 
ages due to thermal fluctuations. We thus introduce a radial 
mask to remove surface fluctuations from the detection pro- 
cess by determining the maximum penetration depth of sur- 
face fluctuation phase singularities for a given temperature; 
this is accomplished by an adaptive routine that is applied to 
the c-field data for each trajectory. 

At short times, the surface excitations are clearly distin- 
guishable from the central vortex. Our numerical algorithm 
for each trajectory is as follows: (a) at each time step we de- 
tect all vortices with radius < Ry. (b) we discard any that 
are further than a small numerical tolerance from the loca- 
tion of the decaying vortex at the previous time step. is 
then adjusted inwards to the smallest radius of any spurious 
vortex signal. The mask radius moves rapidly inward until it 
eliminates all surface fluctuations from the vortex detection 
region. We take this to be the largest vortex core displace- 
ment radius at which the vortex is stifl resolvable, Ry, and 
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FIG. 5. Column densities and phase slices showing SPGPE vortex evolution in a system with aspect ratio of 4:1 at temperatures (a) T/Tc = 
0.78, and (b) T/Tc = 0.93. The vortex is most easily identified from the phase singularity (circle), and the center of the condensate is marked 
by a cross (+). Other parameters are for a gas of 5 x 10^ *^Rb atoms in a harmonic trap with (ujr,uJz) = 27r(12.4, 49.7)Hz. 



use it to define the vortex lifetime in a single SPGPE trajec- 
tory ty = ma,x{t : r.y < i?„).The specific value of i?„ is 
temperature, and, less directly, number dependent. At zero 
temperature Ry ~ Rtf', for T — > Tc, Ry — > 0, demonstrat- 
ing the shrinking phase coherent region inside the BEC. As 
the SPGPE is a stochastic theory, in what follows we compute 
the mean vortex lifetime, t = (ty), over a set of trajectories 
for each temperature considered. For simplicity, we hereafter 
refer to the trajectory-averaged t as the vortex lifetime. We 
find that this procedure gives a time for vortex death that cor- 
responds to a specific value of angular momentum per particle 
0.3fi) that is independent of the system temperature. 



III. RESULTS 

To understand the different decay mechanisms we first 
compare the dynamics of a centrally located vortex by study- 
ing single trajectories within PGPE, DPGPE, and SPGPE the- 
ories. 



A. Comparison of models of vortex decay 



Using the initial finite T central vortex states (see section 



II B I, we compare subsequent evolution with respect to the 



three approaches. For the PGPE case, we neglect all inter- 
actions with the above cutoff region in the SPGPE. For the 
DPGPE case we neglect only the noise, and we also introduce 
a small initial offset of the vortex. The offset is necessary 
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T/T, 


Theory 


t (sec) 1:1 


t (sec) 4:1 




PGPE 


oo 


oo 


0.78 


DPGPE 


33.7 


69.4 




SPGPE 


10.7 


20.4 




PGPE 


2.6 


3.5 


0.93 


DPGPE 


21.0 


53.9 




SPGPE 


1.5 


2.0 



TABLE II. Comparison of vortex lifetime in different theories for a 
vortex in a Bose gas with spherical (1:1) and oblate (4:1) geometries. 



25 r 



20- 



15- 



10- 



5- 



O Lj^/oJr = 4 
iJ^/u!,. = 1 



because a central vortex is metastable, and hence evolution 
according to the purely damped PGPE will not exhibit vor- 
tex decay unless there is sufficient symmetry breaking. Initial 
fluctuations and noise in the PGPE and SPGPE approaches 
provide sufficient symmetry breaking to cause decay, but the 
initial fluctuations in the DPGPE treatment are insufficient. 
For the DPGPE and SPGPE, cases the value of 7 is deter- 



0.76 0.78 0.8 0.82 0.84 0.86 O.f 



0.9 0.92 0.94 



mined by Eq. ( 14 1. 

A comparison of individual trajectories for the DPGPE and 
PGPE approaches is presented in Figure |4] for a system in a 
4:1 trapping geometry with (ujrjUJz) = 27r x (12.4, 49.7)Hz. 
In Fig.|4](a) and (b) the initial temperature is T = 34.0 nK. 
In the DPGPE approach shown in Fig. ^a), initial thermal 
fluctuations are quickly damped out, the field density be- 
comes smooth, and the vortex then decays according to fa- 
miliar damped GP evolution. The small value of hj/ksT = 
7.1 X 10^^ leads to a long lifetime of t — 69.4s. Fig- 
ure[4];b) shows PGPE evolution for the same system. The vor- 
tex evolves into a quasi-equilibrium state of precession about 
the z-axis at an approximately fixed radius. At this value of 
temperature, the angular momentum in the initial vortex con- 
figuration is too large to be completely transferred to the non- 
condensate fraction of the classical field, inhibiting vortex de- 
cay. 

In Fig. |4|c) the initial state has temperature T — 64.9 nK 
and thus contains a significantly larger non-condensate frac- 
tion. After a short time (3.5s), the vortex irreversibly leaves 
the condensate through exchange of angular momentum with 
non-condsate degrees of freedom. Figure |5j a) shows SPGPE 
trajectory results for the same parameters as Fig. Qa), (b). 
While the vortex did not decay in PGPE and decayed slowly 
in DPGPE, the combined influence of classical fluctuations in 
the c-field and reservoir interactions leads to a much shorter 
lifetime of ~ 20s in the 4: 1 geometry. Figure |5|b) shows an 
SPGPE trajectory for the parameters of Fig.gc). The SPGPE 
lifetime of 2s is significantly shorter than the 3.5s lifetime of 
the PGPE, which is shorter than the DPGPE lifetime of 54s. 
The results for lifetimes are summarized in Table|ll) where we 
also show lifetimes for a spherical trap. We note that the tem- 
perature dependence of the lifetime of a vortex evolved under 
the DPGPE is much weaker than under SPGPE evolution; the 
SPGPE lifetime changes by an order of magnitude between 
the two temperatures considered. 

While we have not yet performed a systematic study of the 



FIG. 6. Lifetime of a single vortex, i, as a function of temperature 
for a system containing 5 x 10^ *^Rb atoms. The trap has fixed 
= 27r X 19.7 so that Tc = 69.7 nK is the same for both 
geometries. 



influence of trap geometry, we note that the main effect of re- 
ducing the oblateness of the trap (note that is held constant) 
appears to be to shorten the lifetime further. Presumably this 
is associated with the accessibility of vortex bending modes or 
Kelvin waves [34J which provide a further dissipative mecha- 
nism that is not available in highly oblate configurations. 



B. SPGPE simulations of vortex decay 

We perform simulations of the SPGPE for a range of tem- 
peratures and calculate the vortex lifetime t. All quantities of 
interest are averaged over an ensemble of trajectories which is 
chosen sufficiently large to give an acceptable level of stochas- 
tic convergence. For each temperature we choose parameters 
in order to keep iVr = 5 x 10^ constant, using the estima- 
tion procedure of Sec. [ll] The resulting /i is used to evaluate 
Eq. ( [T4| to determine the strength of dissipation and noise in 
the SPGPE. 

The mean vortex lifetime for 4:1 and 1:1 geometries is 
shown in Fig. [6] For both geometries the geometric mean and 
temperature are iIi = 27r x 19.7 s^^, and Tc = 69.7 nK 
respectively. Rather surprisingly, we observe that the vor- 
tex lifetime is linear in T/Tc at fixed Nt. This result may 
not be totally unexpected since the leading dependence of the 
damping on temperature is given by the prefactor of the di- 
mensionless damping rate h'y/kBT ^ Ama?kBT /ttH^ . How- 
ever, for constant A^t the temperature variation leads to large 
changes in chemical potential, and thus nonlinear changes in 
7 L12J. The characteristic size of the condensate also shrinks 
with increasing T [decreasing /z, see Fig.|3](a)], so the vortex 
is not only more mobile, it has a shorter distance to traverse 
(e.g. the condensate radius in Fig. |5ja) is ^ 15/im, while 
in Fig. [5](b) it has reduced to ^ 10/im). We note that the 
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slope of t increases with iOzjiOr- We again associate this in- 
creased dissipation with the accessibihty of vortex bending 
modes for the more spherical system. This connection can be 
made because we have held all parameters of the system con- 
stant, apart from temperature and the ratio LOz/oJr- Spurious 
factors that would change the vortex decay in an uncontrolled 
way, such as changing A'x or oj, have been eliminated. 

Figure |7] shows the angular momentum per particle of the 
C region during vortex decay. To identify regimes of dif- 
ferent qualitative behavior we graph the data in units of the 
vortex lifetime t calculated for each temperature. The point 
t = i nearly coincides with the intersection point of all curves, 
namely where {Lz) /Nq « 0.3. This corresponds to the angu- 
lar momentum remaining when the vortex is indistinguishable 
from thermal fluctuations at the condensate boundary. We also 
plot {Lz) /Nq for the overdamped limit of the DPGPE, which 
drops to zero att — t; this corresponds to the time when the 
vortex reaches the Thomas-Fermi radius of the BEC. Compar- 
ing with the overdamped limit, we see that for low tempera- 
tures, the form of the decay in the SPGPE approaches that of 
purely dissipative decay governed by the DPGPE. The main 
effect of the additional thermal fluctuations in the SPGPE is to 
cause rapid loss of angular momentum at short times, corre- 
sponding to destabilizing the vortex from its initial position 
at the origin. Thus, thermal fluctuations render the vortex 
metastable. 

As temperature increases, the angular momentum becomes 
increasingly persistent for t < t. The angular momentum 
develops a broad tail at long times associated with a broader 
distribution of vortex radii. We refer to this high temperature 
behavior as the regime of diffusive decay, where the physics of 
vortex decay is dominated by the diffusive influence of ther- 
mal fluctuations. 



IV. CONCLUSIONS 

We have presented a scheme to conveniently estimate con- 
sistent parameters in the regime of validity of the SPGPE for 
systems near thermal equilibrium. The chemical potential is 
calculated using a Hartree-Fock calculation of the C-region 
occupation. The cutoff energy e^ut is an essential feature of 
the theory used to consistently separate the coherent and inco- 
herent regions of the system. It is determined by finding the 
single particle cutoff that bounds the same number of quantum 
states as contained in the Hartree-Fock treatment of C, and 
implemented numerically using the projector in the SPGPE. 
This corrects the single particle cutoff for the compression of 
energy levels caused by the presence of the condensate. By 
finding consistent SPGPE parameters, we are able to set the 
value of the damping constant 7 through Eq. 
determine all physical parameters of the theory. 



14 1, and thus 




Comparing three c-field theories capable of describing vor- 
tex decay in different levels of approximation, we found that 
the SPGPE theory leads to the fastest decay in all cases con- 
sidered. While this is expected as the SPGPE contains the 
other theories as special cases, in some cases the differences 
between theories are quite dramatic. At the lower temperature 



FIG. 7. The z component of the averaged angular momentum per 
particle as a function of t/t, where t is the vortex lifetime, for a 
system with an aspect ratio of 4: 1 as found from an ensemble average 
of SPGPE trajectories. The case of purely dissipative decay is given 
by the dotted line. Note the wide contrast in the form of the vortex 
decay where for T/Tc = 0.78 we see dissipation dominating the 
decay, while for higher temperatures diffusive decay becomes more 
prevalent. 



considered, complete decay of the vortex was not observed 
in the PGPE theory in either the 4:1 or 1:1 geometries. Vor- 
tex decay was observed at the higher temperature, providing 
another example of ergodicity emerging from purely Hamil- 
tonian evolution. For all cases where a finite lifetime is ob- 
served, the DPGPE led to slower decay that the PGPE. This 
is caused by the rapid damping in DPGPE of thermally occu- 
pied classical degrees of freedom in the pure PGPE theory. In 
order for the damping to give sensible results, it is thus vital 
to include the associated stochastic element of the reservoir 
interaction, leading to the SPGPE theory. When the temper- 
ature is high enough to give a finite lifetime in PGPE theory, 
the calculated value is within a factor of two of the predictions 
of the SPGPE. 

Applying the SPGPE theory to vortex decay over a range 
of temperatures we observe that the vortex lifetime is lin- 
early dependent on T/Tc at fixed iVr, and varies by approx- 
imately an order of magnitude over the temperature range 
O.TSTc < T < 0.93T;. The shortest lifetimes, calculated for 
5x 10^ «^Rb atoms at T/Tc = 0.93 in 4:1 and 1:1 geometries, 
are i — 2s and fc=1.5s respectively. At the high temperatures 
corresponding to these short lifetimes, the angular momentum 
gives a clear indication of diffusive decay of the vortex, evi- 
denced by a long tail for t < t. 

We note that our results for the vortex lifetime also pro- 
vide a lower bound for the lifetime of a persistent current of 
unit winding number, as any additional vortex pinning po- 
tential will further extend the vortex lifetime. Our results 
have been calculated using a theory that we expect to be a 
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very good description of the high temperature Bose gas in 
g'Masi-equilibrium, without any fitting parameters. We are 
thus optimistic that our predictions will be of quantitative 
value in comparing with future experimental measurements 
of the lifetimes of vortices and persistent currents. Future 
work will focus on the role of the scattering term in the full 
SPGPE f7' T^ '23l, and a more detailed treatment of the inco- 
herent region dynamics. 

We are grateful to T. M. Wright for stimulating discussions 
regarding this work. We acknowledge support from the New 
Zealand Foundation for Research, Science and Technology 
contracts NERF-UOOX0703, and UOOX0801. 



Appendix A: Hartree-Fock estimation of SPGPE parameters 

In this appendix we give details of our method for deter- 
mining consistent control parameters for the SPGPE theory /i, 
Ecut^ as a function of temperature T, for a system constrained 
to have fixed total particle number iVr. The method is based 
on a Hartree-Fock treatment of the trapped Bose gas. 

In the Thomas-Fermi regime, the condensate density is 
given by 

I /i— K(x)<OJ 

where /i sets the condensate occupation Nq — J nQ{x)(Px. 
Using Eq. (21 1 to describe thermal excitations, the Hartree- 
Fock density of states for the system 



PHF(e) 



-6{e- EiiF{K,p)) . 



(A2) 



can be used to determine the total atom number for a given 
/i. For the harmonic trap, this can be evaluated to give the 
analytical expression of the form 



PHF(e) = ^ [/-(e) +/+(£)] 

where /-(e) and /+(e) are defined by 

^-i^)={—^ ^ 



(A3) 



-log(a; + w_) 



x—^pi^rfhw 

x—^ max{0,a_ } 



(A4) 



/+(e) 



-sm 




a;— ^a-f 



(A5) 



and where a± = 2(e ± /ixp)/^, and u± — y^a±^ x"^. 



The total number of atoms in the system is given by 



where 



dt PHF(e)n-BE(e) + Nq 



nBE(e) 



1 



(A6) 



(A7) 



By specifying the total number of atoms in the system equa- 



tion ( A6 1 can be solved numerically to give A^o^ and hence /i, 
using the Thomas-Fermi chemical potential for the harmonic 
trap 



2/5 



(A8) 



where a = ^hlrriib. For temperatures above Tc, /OHF(e) 
becomes the standard single particle density of states. The 
SPGPE parameter /i(r) can then be found easily for temper- 



atures above the transition using ( A6 1, with iVo = 0. 

We also need to determine ecut so that the highest energy 
C-region modes satisfy our occupation condition ricut- How- 
ever, some care is required to identify the value of e^ut for 
use in the SPGPE. The SPGPE simulations are most conve- 
niently performed using a spectral numerical tecnique based 
on the harmonic oscillator eigenstates of the trap |35|. Thus 
Ecut is chosen as the highest oscillator eigenstate included in 
the C region. The actual physical energy of the corresponding 
oscillator state in the SPGPE c-field is then shifted by inter- 
actions. The simplest way to connect ecut with CcutHF is to 
require that the number of modes in the C-region is the same 
in either representation. 

The Hatree-Fock density of states takes energies relative to 
the condensate, so from ( |A7[ ) with /i — we have 



CcutHF — fce/'ln 1 



1 

ncut 



(A9) 



valid for T < where CcutHF is the energy set by the HE 
density of states at mean occupation ricut- The relevant num- 
ber of single particle states is given by integrating phf up to 
EcutHF since the condensate compresses the low lying levels 
together without changing their number, ecut is then found by 
requiring that 



rfe PHF(e) 



&PHo(e), (AlO) 



where eq — ti{2uJr + i^z)/"^, and paoi^) — e^/2(/ia;)^. For 
temperatures above the transition, the single particle density 
of states is applicable so ecut can be found directly from the 
Bose-Einstein distribution 



CcutHF = Ecut = kBTln 1 



1 

ncut 



H. (All) 



This set of equations can be rapidly evaluated to give fJ.{T), 
and ecut{T) [see Fig. [3). 
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